!! Copyright (C) 2010,2011,2012  Luca Bonaventura, MOX, Polimi
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> \brief
!!
!! Splitting methods with exponential time integration
!!
!! \n
!!
!! This module provides additional integrators with respect to those
!! implemented in \c mod_expo, which include a splitting of the
!! right-hand side. The main reference for these methods is <a
!! href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo, Saad,
!! J. Sci. Comp., 1998]</a>.
!!
!! We discuss here some aspects which are relevant for all the
!! considered methods. Let us thus consider
!! \f{displaymath}{
!!  y' = f(t,y) = f_1(t,y) + f_2(t,y)
!! \f}
!! where \f$f_1\f$ will be treated with exponential integration and
!! \f$f_2\f$ will be treated explicitly. Following what is done in \c
!! mod_expo, let us linearize the first term
!! \f{displaymath}{
!!  f_1(t,y) = f_1(t^n,y^n) + B(y-y^n) + R(t,y),
!! \f}
!! with \f$\displaystyle B=\frac{df_1}{dy}(t^n,y^n)\f$, and consider
!! the problem
!! \f{displaymath}{
!!  \left\{
!!  \begin{array}{rcll}
!!  y' -By & = & g^n + \tilde{f}_2(t,y), & t\in[t^n\,,t] \\
!!  y(t^n) & = & y^n
!!  \end{array}
!!  \right.
!! \f}
!! with \f$g^n = f_1(t^n,y^n)-By^n\f$ and \f$\tilde{f}_2(t,y)
!! =f_2(t,y)+R(t,y)\f$. The exact solution to this equation is
!! \f{displaymath}{
!!  y(t) = e^{B(t-t^n)}\left(
!!   y^n + \int_{t^n}^{t}e^{-B(s-t^n)}\left( g^n + \tilde{f}_2(s,y(s))
!!   \right)ds \right),
!! \f}
!! so that we have the (formal) time step 
!! \f{displaymath}{
!!  y^{n+1} = y^n + 
!!   \int_{t^n}^{t^{n+1}}e^{-B(s-t^{n+1})}\left( f_1(t^n,y^n) +
!!   \tilde{f}_2(s,y(s)) \right)ds,
!! \f}
!! where we have used
!! \f{displaymath}{
!!  e^{\Delta t\,B} = I + B\int_0^{\Delta t}e^{-B(s-\Delta t)}ds
!!  = I + B\int_{t^n}^{t^{n+1}}e^{-B(s-t^{n+1})}ds.
!! \f}
!! \note It is often possible to assume \f$\tilde{f}_2\approx f_2\f$,
!! either because \f$f_1\f$ is linear or because the nonlinearity can
!! be neglected <em>within the time step</em>. Formally however,
!! neglecting the nonlinearity in \f$f_1\f$ limits the order of the
!! method to second order.
!!
!! \section PCM Predictor corrector methods
!!
!! Predictor-corrector methods are discussed in <a
!! href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo, Saad,
!! J. Sci. Comp., 1998]</a>, and the main idea is to approximate with
!! the predictor step the value \f$\tilde{f}_2(s,y(s))\f$ appearing in
!! the formal time step, in order to be able to compute the integral
!! (corrector step). In fact, assuming \f$\tilde{f}_2(s,y(s))\approx K
!! = \tilde{f}_2(\bar{t},\bar{y})\f$, for some fixed time level
!! \f$(\bar{t},\bar{y})\f$ and evaluating the integral yields
!! \f{displaymath}{
!!  y^{n+1} = y^n + 
!!   \Delta t\,\varphi(\Delta t\,B)\left(f_1(t^n,y^n)+K \right).
!! \f}
!!
!! The simplest algorithms of this family are obtained using
!! \f$\tilde{f}_2(s,y(s))\approx f(t^n,y^n)\f$ (this method does not
!! require a predictor step) and \f$\tilde{f}_2(s,y(s))\approx
!! f(t^n+\frac{\Delta t}{2},y^{n+\frac{1}{2}})\f$. The second method
!! corresponds to Algorithm 2.1 of <a
!! href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo, Saad,
!! J. Sci. Comp., 1998]</a>. The advantage of these two methods is
!! that \f$\tilde{f}_2(s,y(s))\f$ is approximated by a constant, thus
!! simplifying the evaluation of the integral.
!!
!! To obtain higher order methods, however, one has to consider (for
!! instance) a polynomial representation of \f$\tilde{f}_2(s,y(s))\f$,
!! so that one has to deal with integrals of the form
!! \f{displaymath}{
!!  \int_{t^n}^{t^{n+1}}s^ke^{-B(s-t^{n+1})}ds
!! \f}
!! for values \f$k\geq0\f$. These integrals can be expressed in terms
!! of the functions \f$\varphi_k\f$ introduced in <a
!! href="http://dx.doi.org/10.1137/0729014">[Saad, SIAM J.  Numer.
!! Anal., 1992]</a>, section 5.1, <a
!! href="http://dx.doi.org/10.1007/s00607-007-0227-1">[Caliari,
!! Computing, 2007]</a>, section 3.1, <a
!! href="http://dx.doi.org/10.1016/j.apnum.2008.03.021">[Caliari,
!! Ostermann, Appl. Num. Math., 2009]</a>, section 2 and <a
!! href="http://dx.doi.org/10.1016/j.procs.2010.04.026">[Tokman,
!! Loffeld, Proc. Comp. Sci., 2012]</a>, section 2.
!!
!! \section GRKM Generalized Runge-Kutta methods
!!
!! In <a href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo,
!! Saad, J. Sci. Comp., 1998]</a> the idea of <em>generalized
!! Runge-Kutta methods</em> proposed in <a
!! href="http://dx.doi.org/10.1137/0704033">[Lawson, SIAM J. Sci.
!! Comp., 1967]</a> is considered, which is based on rewriting the ODE
!! system in a more tractable way at the continuous level and then
!! applying a classical Runge-Kutta method to this modified (scaled)
!! problem. In practice, one introduces the scaled variable
!! \f$z(t)=e^{-B(t-t^n)}y(t)\f$ and rewrites the problem as
!! \f{displaymath}{
!!  \left\{
!!  \begin{array}{rcll}
!!  z' & = & \displaystyle
!!    \underbrace{e^{-B(t-t^n)}\left( g^n +
!!    \tilde{f}_2(t,e^{B(t-t^n)}z) \right)}_{g(t,z)},& t\in[t^n\,,t]
!!    \\
!!  z(t^n) & = & y^n.
!!  \end{array}
!!  \right.
!! \f}
!! Applying a standard RK method to such a problem results in an
!! exponential RK formulation for the original equation.
!!
!! In terms of the scaled variable this amounts to stages
!! \f{displaymath}{
!!  Z_i = z^n + \Delta t\sum_j a_{ij}K_j
!! \f}
!! and update
!! \f{displaymath}{
!!  z^{n+1} = z^n + \Delta t\sum_j b_jK_j,
!! \f}
!! with
!! \f{displaymath}{
!!  K_j = g(t^n+c_j\Delta t,Z_j).
!! \f}
!! In terms of the original variable, this becomes
!! \f{displaymath}{
!!  Y_i = e^{c_i\Delta t\,B} \left(
!!  y^n + \Delta t\sum_ja_{ij}e^{-c_j\Delta t\,B}(
!!  g^n + \tilde{f}_2(t^n+c_j\Delta t,Y_j) )
!!  \right)
!! \f}
!! and
!! \f{displaymath}{
!!  y^{n+1} = e^{\Delta t\,B} \left(
!!  y^n + \Delta t\sum_jb_{j}e^{-c_j\Delta t\,B}(
!!  g^n + \tilde{f}_2(t^n+c_j\Delta t,Y_j) )
!!  \right).
!! \f}
!! This formalism has the inconvenience that the RK stages are given
!! in terms of the matrix exponential, rather than in terms of
!! \f$\varphi(\Delta t\,B)\f$, so that their implementation on top of
!! the exponential Euler method is not straightforward.
!!
!! The simplest method of this family is obtained using the explicit
!! Euler method to integrate the scaled equation, which provides
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!  y^{n+1} & = & \displaystyle e^{\Delta t\,B} \left(
!!  y^n - \Delta t\,By^n + \Delta t\,f(t^n,y^n)
!!  \right) \\[3mm]
!!  & = & \displaystyle
!!  \underbrace{y^n+\Delta t\,\varphi(\Delta t\,B)f(t^n,y^n)}_{ {\rm
!!  exp.\,\, Rosenbrock\,\, 2nd\,\, order}} + \underbrace{\Delta
!!  t\left(\varphi(\Delta t\,B)-e^{\Delta t\,B}\right)}_{
!!  \mathcal{O}(\Delta t^2)} \left( By^n-f(t^n,y^n)\right).
!!  \end{array}
!! \f}
!! This shows that even in the case \f$f_2=0\f$, where \f$B\f$ is the
!! Jacobian of \f$f\f$, the method does not coincide with the lowest
!! order exponential Rosenbrock method discussed in \c mod_expo
!! (unless \f$f\f$ is linear, in which case both methods yield the
!! exact solution).
!!
!! The second method of this family is obtained using the Heun method
!! for the scaled equation, and is described as Algorithm 2.2 in
!! <a href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo,
!! Saad, J. Sci. Comp., 1998]</a>. The final update can be written as 
!! \f{displaymath}{
!!  y^{n+1} = 
!!  \underbrace{y^n+\Delta t\,\varphi(\Delta t\,B)f(t^n,y^n)}_{ {\rm
!!  exp.\,\, Rosenbrock\,\, 2nd\,\, order}} + \underbrace{\Delta
!!  t\left(\varphi(\Delta t\,B)-\frac{e^{\Delta t\,B}+1}{2}\right)}_{
!!  \mathcal{O}(\Delta t^3)} \left( By^n-f(t^n,y^n)\right) + X,
!! \f}
!! where \f$X = \frac{\Delta t}{2}\left( \tilde{f}_2(t^n+\Delta t,Y_2)
!! - f_2(t^n,y^n)\right)\f$ and \f$Y_2\f$ has the same expression
!! obtained for the final update of the explicit Euler method. Again,
!! it can be seen that even for \f$f_2=0\f$ this is indeed a new
!! method, different from the exponential Rosenbrock one.
!!
!! \subsection GRKMphi An alternative formulation
!!
!! To obtain a similar method in terms of the function\f$\varphi\f$ we
!! can proceed as follows. Let us consider the affine problem
!! \f{displaymath}{
!!  y' = By+q,
!! \f}
!! for a constant \f$q\f$, for which we would like to recover the
!! exact solution
!! \f{displaymath}{
!!  y(t) = y^n + B^{-1}\left( e^{B(t-t^n)}-I\right)\left(
!!  By^n+q\right).
!! \f}
!! Two possible choices for the scaled variable are now
!! \f{displaymath}{
!!  z(t) = B\left(e^{B(t-t^n)}-I\right)^{-1}\left(y(t)-y^n\right)
!!  \approx By^n+q,
!! \f}
!! and
!! \f{displaymath}{
!!  z(t) = \underbrace{
!!  (t-t^n)B\left(e^{B(t-t^n)}-I\right)^{-1}}_{
!!  \left(\varphi(B(t-t^n))\right)^{-1}}
!!  \left(y(t)-y^n\right) \approx
!!  (t-t^n)\left(By^n+q\right),
!! \f}
!! since any Runge-Kutta method is exact for constant or linear
!! solutions (the symbol \f$\approx\f$ must be understood as \f$=\f$
!! when \f$y\f$ is the solution of the affine problem). Taking the
!! first choice, we have
!! \f{displaymath}{
!!  z'(t) = \frac{1}{t-t^n}\varphi^{-1}\left[ f(t,\varphi\,z+y^n)- 
!!  e^{B(t-t^n)} z\right],
!! \f}
!! where the argument \f$B(t-t^n)\f$ of \f$\varphi\f$ is omitted for
!! readability, while the second choice yields
!! \f{displaymath}{
!!  z'(t) = \varphi^{-1}\left[ f(t,\varphi\,z+y^n)- \frac{1}{t-t^n}\left(
!!  e^{B(t-t^n)}-\varphi\right) z\right].
!! \f}
!! In both cases, discretizing the scaled problem with the explicit
!! Euler method results in the exponential Rosenbrock integrator.
!!
!! \note The difference with the generalized RK methods of 
!! <a href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo,
!! Saad, J. Sci. Comp., 1998]</a> is that those are exact for linear
!! problems \f$y'=By\f$, while this alternative form is exact for
!! affine problems \f$y'=By+q\f$. However, one must notice that for
!! higher order RK method this alternative form can become rather
!! impractical.
!!
!! \section ERBW Exponential Rosenbrock-Wanner schemes
!!
!! Exponential Rosenbrock methods have been introduced in <a
!! href="http://dx.doi.org/10.1137/S0036142995280572">[Hochbruck,
!! Lubich, SIAM J. Numer. Anal., 1997]</a> and <a
!! href="http://dx.doi.org/10.1137/S1064827595295337">[Hochbruck,
!! Lubich, Selhofer, SIAM J. Sci. Comp., 1998]</a>. The main idea
!! behind these methods is the following:
!! <ul>
!!  <li> starting from a diagonally implicit RK method, linearize the
!!  implicit problem introducing the Jacobian of the right-hand-side,
!!  and to obtain a more general method allow for linear combinations
!!  of the Jacobian terms with coefficients \f$\gamma_{ij}\f$; the
!!  resulting method is a Rosenbrock method, described for instance in
!!  <a href="http://dx.doi.org/10.1007/978-3-642-05221-7">[Hairer,
!!  Wanner, Solving Ordinary Differential Equations II]</a>, Chapter
!!  IV.7;
!!  <li> generalize the method by considering an arbitrary
!!  approximation of the Jacobian matrix; the resulting method is a
!!  Rosenbrock-Wanner scheme, described again in <a
!!  href="http://dx.doi.org/10.1007/978-3-642-05221-7">[Hairer,
!!  Wanner, Solving Ordinary Differential Equations II]</a>, Chapter
!!  IV.7;
!!  <li> finally, as discussed in <a
!!  href="http://dx.doi.org/10.1137/S0036142995280572">[Hochbruck,
!!  Lubich, SIAM J. Numer. Anal., 1997]</a>, Section 5, replace the
!!  matrix \f$(I-\gamma\Delta t\,B)^{-1}\f$ of the linear system to be
!!  solved at each stage with the exponential matrix
!!  \f$\varphi(\gamma\Delta t\,B)\f$ (and possibly recompute the
!!  coefficients to ensure that the order conditions are verified, see
!!  section 3.1 and 3.2 of <a
!!  href="http://dx.doi.org/10.1137/S1064827595295337">[Hochbruck,
!!  Lubich, Selhofer, SIAM J. Sci. Comp., 1998]</a>).
!! </ul>
!! \note The substitution \f$(I-\gamma\Delta t\,B)^{-1}\mapsto
!! \varphi(\gamma\Delta t\,B)\f$ is equivalent to the assumption that
!! during the \f$i\f$-th RK step the solution evolves according to
!! \f$y'=By + q_i\f$, with initial condition \f$y^n\f$, for some
!! constant \f$q_i\f$.
!!
!! These methods have the convenient feature that they can be
!! implemented in terms of the sole matrix function \f$\varphi\f$,
!! which simplifies implementing them on top of the methods described
!! in \c mod_expo. The general form is (neglecting the partial
!! derivative with respect to \f$t\f$ in case of nonautonomous
!! systems)
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  Y_i & = & \displaystyle y^n +\Delta t\,
!!  \sum_j\alpha_{ij}K_j \\[4mm]
!!  y^{n+1} & = & \displaystyle y^n +\Delta t\,
!!  \sum_jb_{j}K_j \\[4mm]
!!  K_i & = & \displaystyle \varphi(\gamma\Delta t\,B)\left(
!!  f(t^n+c_i\Delta t,Y_i) + \Delta t\,B\sum_j\gamma_{ij}K_j\right).
!!  \end{array}
!! \f}
!! The simplest method of this family, corresponding to the explicit
!! Euler method, coincides with the method implemented in \c mod_expo.
!!
!! \section others Further references
!!
!! We collect in this section some references to related methods.
!! <ul>
!!  <li> <a href="http://dx.doi.org/10.3934/dcds.2014.34.1079">[Lee,
!!  Engquist]</a>, method for multiscale (in time) ODEs
!!  <li> <a href="http://dx.doi.org/10.1007/s10915-013-9718-8">[Misun
!!  Min, Paul Fischer]</a> exponential integrators for DG
!!  <li> <a href="http://dx.doi.org/10.1016/j.cam.2012.09.038">[ J.
!!  Loffeld, M. Tokman]</a> comparison implicit vs. exponential - the
!!  exponential seems more efficient in various cases
!!  <li> <a href="http://dx.doi.org/10.1137/110820191">[Botchev, M.,
!!  Grimm, V., and Hochbruck, M.]</a> restarting the Krylov method for
!!  the computation of the matrix exponential.
!!  <li> <a href="http://dx.doi.org/10.1137/080717717">[Hochbruck,
!!  M., Ostermann, A., and Schweitzer, J., 2009]</a> overview for
!!  exponential Rosenbrock methods.
!! </ul>
module mod_spexpo

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp, wp_p, wp_r

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_logical, mpi_integer, wp_mpi, &
   mpi_comm_world, &
   mpi_bcast

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_time_integrators_base, only: &
   mod_time_integrators_base_initialized, &
   c_ode, c_ods, c_tint, t_ti_init_data, t_ti_step_diag

 use mod_expo, only: &
   mod_expo_initialized, &
   t_erb2kry, t_erb2lej, &
   c_expo, erb2kry_step_internals, erb2lej_step_internals

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_spexpo_constructor, &
   mod_spexpo_destructor,  &
   mod_spexpo_initialized, &
   t_pcexpo

 private

!-----------------------------------------------------------------------

! Module variables

 ! public members
 logical, protected ::               &
   mod_spexpo_initialized = .false.

 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_spexpo'

!-----------------------------------------------------------------------

 !> Data specific to the split exponential methods
 !!
 !! The main purpose of these fields is to allow calling the functions
 !! defined in \c mod_expo, providing them some storage to allocate
 !! the required working arrays and data.
 type t_spexpo_data
  !> erb2 integrator to call <tt>erb2*_step_internals</tt>
  class(c_expo), allocatable :: erb2
  !> Object used to recover the diagnostics from
  !! <tt>erb2*_step_internals</tt>
  type(t_ti_step_diag) :: erb2_step_diag
  !> Pointer to the selected <tt>erb2*_step_internals</tt>
  procedure(erb2kry_step_internals), pointer, nopass :: erb2step=>null()
  !> MPI specific data
  integer :: mpi_id = 0
  integer :: mpi_comm = mpi_comm_world
 end type t_spexpo_data

 !> Predictor-corrector exponential, second order
 !!
 !! This is the PC-EXPO method described in algorithm 2.1 in <a
 !! href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo, Saad,
 !! J. Sci. Comp., 1998]</a>, implemented following the framework
 !! described in the \ref PCM "module documentation".
 !!
 !! For this method, the integral defining the time increment is
 !! evaluated with a predictor-corrector method. Namely:
 !! <ul>
 !!  <li> predictor: compute an intermediate time level assuming
 !!  \f$f_2(s,y(s))\approx f_2(t^n,y^n)\f$
 !!   \f{displaymath}{
 !!    \begin{array}{rcl}
 !!    \displaystyle
 !!    y^{n+\frac{1}{2}} & = & \displaystyle y^n + 
 !!   \int_{t^n}^{t^n+\frac{\Delta t}{2}} e^{-B(s-t^n-\frac{\Delta
 !!   t}{2})} \left( f_1(t^n,y^n) + f_2(t^n,y^n) \right)ds \\[3mm]
 !!    & = & \displaystyle
 !!    y^n + \frac{\Delta t}{2}\varphi\left(\frac{\Delta
 !!    t}{2}B\right)\,f(t^n,y^n)
 !!    \end{array}
 !!   \f}
 !!  <li> corrector: compute the whole time step with the
 !!  approximation
 !!  \f$f_2(s,y(s))\approx f_2(t^n+\frac{\Delta
 !!  t}{2},y^{n+\frac{1}{2}})\f$
 !!   \f{displaymath}{
 !!    \begin{array}{rcl}
 !!    \displaystyle
 !!    y^{n+1} & = & \displaystyle y^n + 
 !!   \int_{t^n}^{t^{n+1}} e^{-B(s-t^{n+1})} \left( f_1(t^n,y^n) +
 !!   f_2\left(t^n+\frac{\Delta t}{2},y^{n+\frac{1}{2}}\right)
 !!   \right)ds \\[3mm]
 !!    & = & \displaystyle
 !!    y^n + \Delta t\,\varphi\left(\Delta
 !!    t\,B\right)\,\left(f_1(t^n,y^n) +
 !!   f_2\left(t^n+\frac{\Delta t}{2},y^{n+\frac{1}{2}}\right)
 !!    \right).
 !!    \end{array}
 !!   \f}
 !! </ul>
 !! Remarks:
 !! <ul>
 !!  <li> both steps correspond to a standard exponential integrator
 !!  where the matrix exponential is computed using the reduced matrix
 !!  \f$B\f$ and a complete right-hand side \f$f\f$
 !!  <li> the two steps use time increment \f$\frac{\Delta t}{2}\f$
 !!  and \f$\Delta t\f$, respectively
 !!  <li> \f$f_2\f$ enters the method defining the vector multiplied
 !!  by the exponential matrix, \f$f_1\f$ defining both such vector
 !!  and the exponential matrix
 !!  <li> the fact that \f$f_2\f$ is linear or not does not affect the
 !!  problem in any sense
 !!  <li> if \f$f_1=0\f$ the method reduces to an explicit midpoint
 !!  method which, for linear problem, coincides with the standard
 !!  Heun method (consider that \f$\varphi(0)=I\f$)
 !!  <li> if \f$f_2=0\f$ the first step is useless, while the second
 !!  one reduces to the usual exponential time integrator
 !!  <li> compared to the standard exponential method, a simpler
 !!  matrix exponential is computed, but <em>two</em> matrix
 !!  exponential per time step are required instead of one; one of
 !!  them though uses a halved time step, and thus presumably less
 !!  iterations.
 !!  <li> the accuracy on \f$f_1\f$ comes from the exponential
 !!  Rosenbrock discretization, which includes information about the
 !!  Jacobian, while the accuracy on \f$f_2\f$ comes from the midpoint
 !!  method, and no information is used concerning the derivatives of
 !!  this function.
 !! </ul>
 !!
 !! The detailed implementation works as follows:
 !! <ul>
 !!  <li> <tt>work(1) = </tt> \f$f(t^n,y^n)\f$
 !!  <li> <tt>work(2) = </tt> \f$f_1(t^n,y^n)\f$
 !!  <li> <tt>unew = </tt> \f$y^{n+\frac{1}{2}}\f$
 !!  <li> <tt>work(1) = </tt> \f$f_2\left(t^n+\frac{\Delta t}{2},
 !!  y^{n+\frac{1}{2}}\right)\f$
 !!  <li> <tt>work(1) += work(2)</tt>
 !!  <li> <tt>unew = </tt> \f$y^{n+1}\f$
 !! </ul>
 !!
 !! \section pcexpo_setup Initialization
 !!
 !> The input argument <tt>init_data\%ii1(1)</tt> is used to select
 !! the "single Euler step" method:
 !! <ul>
 !!  <li> 1 \f$\mapsto\f$ \c erb2kry
 !!  <li> 2 \f$\mapsto\f$ \c erb2lej
 !! </ul>
 !! The remaining parameters are passed to the selected single Euler
 !! step method.
 type, extends(c_tint) :: t_pcexpo
  type(t_spexpo_data) :: data
 contains
  procedure, pass(ti) :: init  => pcexpo_init
  procedure, pass(ti) :: step  => pcexpo_step
  procedure, pass(ti) :: clean => pcexpo_clean
 end type t_pcexpo

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_spexpo_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized             .eqv..false.) .or. &
       (mod_kinds_initialized                .eqv..false.) .or. &
       (mod_mpi_utils_initialized            .eqv..false.) .or. &
       (mod_state_vars_initialized           .eqv..false.) .or. &
       (mod_time_integrators_base_initialized.eqv..false.) .or. &
       (mod_expo_initialized                 .eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_spexpo_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_spexpo_initialized = .true.
 end subroutine mod_spexpo_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_spexpo_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_spexpo_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_spexpo_initialized = .false.
 end subroutine mod_spexpo_destructor

!-----------------------------------------------------------------------

 !> Splitted exponential initialization
 subroutine pcexpo_init(ti,ode,dt,t,u0,ods,init_data,step_diag)
  class(c_ode),  intent(in)    :: ode
  real(wp),      intent(in)    :: dt, t
  class(c_stv),  intent(inout) :: u0
  class(c_ods),  intent(inout) :: ods
  class(t_pcexpo), intent(out) :: ti
  type(t_ti_init_data), optional, intent(in)  :: init_data
  type(t_ti_step_diag), optional, intent(out) :: step_diag
 
  integer :: ii1
  character(len=*), parameter :: &
    this_sub_name = 'pcexpo_init'

   ti%dt = dt

   if(present(init_data)) then
     ii1 = init_data%ii1(1)
   else
     ii1 = 1 ! default: Krylov method
   endif

   CaseName: select case(ii1)
    case(1)
     allocate(t_erb2kry::ti%data%erb2)
     ti%data%erb2step => erb2kry_step_internals
    case(2)
     allocate(t_erb2lej::ti%data%erb2)
     ti%data%erb2step => erb2lej_step_internals
    case default
     call error(this_sub_name,this_mod_name,    &
            'Unknown value in "init_data%ii1(1)"')
   end select CaseName
   call ti%data%erb2%init( ode,dt,t,u0,ods,init_data, &
                           ti%data%erb2_step_diag )

   allocate( ti%work(3) , source=u0 )

   ! Diagnostics
   if(present(step_diag)) then
     step_diag%name_d1 = 'ti_conv'
     ! These fields are defined by the single Euler step method
     step_diag%name_d2 = ti%data%erb2_step_diag%name_d1
     step_diag%name_d3 = ti%data%erb2_step_diag%name_d2
     allocate(step_diag%d1(4), step_diag%d2(0,0), step_diag%d3(0,0,0))
   endif

 end subroutine pcexpo_init

!-----------------------------------------------------------------------

 !> Predictor-corrector exponential, second order
 subroutine pcexpo_step(unew,ti,ode,t,unow,ods,step_diag)
  class(t_pcexpo), intent(inout) :: ti
  class(c_ode),  intent(in)    :: ode
  real(wp),      intent(in)    :: t
  class(c_stv),  intent(in)    :: unow
  class(c_ods),  intent(inout) :: ods
  class(c_stv),  intent(inout) :: unew
  type(t_ti_step_diag), optional, intent(inout) :: step_diag

  logical :: master
  integer :: ii(2,3)
  real(wp), allocatable :: tmp_d1(:), tmp_d2(:,:)

   master = ti%data%mpi_id.eq.0

   ! Predictor step
   call ode%rhs(ti%work(1),t,unow,ods,(/2,2/)) ! f2
   call ode%rhs(ti%work(2),t,unow,ods,(/1,2/)) ! f1
   call ti%work(1)%incr(ti%work(2))            ! f = f1+f2
   call ti%data%erb2step(unew,ods ,                      &
     ti%data%erb2%work,ti%data%erb2%data ,               &
     ode,t , unow,ti%work(2),ti%work(1) , ti%dt/2.0_wp , &
     ti%work(3) , ti%data%erb2_step_diag , (/1,2/))
   call unew%incr(unow)

   if(present(step_diag)) then
     ! Store the diagnostics of the first step
     step_diag%d1(1) = l2r(ti%data%erb2_step_diag%max_iter)
     step_diag%d1(2) = real(ti%data%erb2_step_diag%iterations,wp)
     if(master) then
       allocate(tmp_d1(size(ti%data%erb2_step_diag%d1)))
       tmp_d1 = ti%data%erb2_step_diag%d1
       allocate(tmp_d2(size(ti%data%erb2_step_diag%d2,1) , &
                       size(ti%data%erb2_step_diag%d2,2)))
       tmp_d2 = ti%data%erb2_step_diag%d2
     endif
   endif

   ! Corrector
   call ode%rhs(ti%work(1),t+ti%dt/2.0_wp,unew,ods,(/2,2/)) ! f2
   call ti%work(1)%incr(ti%work(2))                         ! f = f1+f2
   call ti%data%erb2step(unew,ods ,                      &
     ti%data%erb2%work,ti%data%erb2%data ,               &
     ode,t , unow,ti%work(2),ti%work(1) , ti%dt        , &
     ti%work(3) , ti%data%erb2_step_diag , (/1,2/))
   call unew%incr(unow)

   if(present(step_diag)) then
     ! Store together the diagnostics of the two steps
     step_diag%d1(3) = l2r(ti%data%erb2_step_diag%max_iter)
     step_diag%d1(4) = real(ti%data%erb2_step_diag%iterations,wp)
     if(master) then
       if(allocated(step_diag%d2)) deallocate(step_diag%d2)
       if(allocated(step_diag%d3)) deallocate(step_diag%d3)
       ii(1,:) = (/ size(tmp_d1) , size(tmp_d2,1) , size(tmp_d2,2) /)
       ii(2,:) = (/ size(ti%data%erb2_step_diag%d1) ,  &
                    size(ti%data%erb2_step_diag%d2,1), &
                    size(ti%data%erb2_step_diag%d2,2) /)
       allocate(step_diag%d2(maxval(ii(:,1)),2))
       step_diag%d2 = -1.0_wp
       step_diag%d2(1:ii(1,1),1) = tmp_d1
       step_diag%d2(1:ii(2,1),2) = ti%data%erb2_step_diag%d1
       allocate(step_diag%d3(maxval(ii(:,2)),maxval(ii(:,3)),2))
       step_diag%d3 = -1.0_wp
       step_diag%d3(1:ii(1,2),1:ii(1,3),1) = tmp_d2
       step_diag%d3(1:ii(2,2),1:ii(2,3),2) = ti%data%erb2_step_diag%d2

       deallocate(tmp_d1,tmp_d2)
     endif
   endif

 contains
  pure function l2r(l)
   logical, intent(in) :: l
   real(wp) :: l2r
    if(l) then
      l2r = 1.0_wp
    else
      l2r = 0.0_wp
    endif
  end function l2r
 end subroutine pcexpo_step

!-----------------------------------------------------------------------

 pure subroutine pcexpo_clean(ti,ode,step_diag)
  class(c_ode),  intent(in)    :: ode
  class(t_pcexpo), intent(inout) :: ti
  type(t_ti_step_diag), optional, intent(inout) :: step_diag

  character(len=*), parameter :: &
    this_sub_name = 'pcexpo_clean'

  call ti%data%erb2%clean( ode , ti%data%erb2_step_diag )
  nullify(ti%data%erb2step)
  deallocate(ti%data%erb2)
  deallocate(ti%work)
  if(present(step_diag)) &
    deallocate(step_diag%d1 , step_diag%d2 , step_diag%d3)

 end subroutine pcexpo_clean

!-----------------------------------------------------------------------

! !> Runge-Kutta exponential, third order
! !!
! !! This is a modified version of the RK3-EXPO method described in
! !! algorithm 2.3 in <a
! !! href="http://dx.doi.org/10.1023/A:1023219016301">[Castillo, Saad,
! !! J. Sci. Comp., 1998]</a>. The modification is introduced to obtain
! !! a method which can be expressed in terms of the \f$\varphi(z)\f$
! !! function, thus simplifying the implementation.
! !!
! !! The main idea is to mimic the RK stages using exponential Euler
! !! steps as discussed in the introduction of this module. The first
! !! stage is simply
! !! \f{displaymath}{
! !!  Y_1 = y^n.
! !! \f}
! !! The second stage is obtained performing an explicit Euler step
! !! with amplitude \f$c_2=\frac{2}{3}\Delta t\f$. Using the
! !! exponential single Euler step, this yields
! !! \f{displaymath}{
! !!  Y_2 = y^n + \frac{2}{3}\Delta t\,\varphi\left( \frac{2}{3}\Delta
! !!  t\, B\right)\left( f_1(t^n,y^n) + \underbrace{f_2(t^n,y^n)}_{K_1}
! !!  \right)
! !! \f}
! !! The third stage uses the trapezoidal rule on the same time
! !! interval as the first stage. This is equivalent to substituting
! !! \f{displaymath}{
! !!  \tilde{f}_2(s,y(s)) \approx 
! !!  \frac{s-t^n}{\Delta t}f_2(t^n,y^n) +
! !!  \frac{t^{n+1}-s}{\Delta t}f_2(t^n+\frac{2}{3}\Delta t,Y_2) =
! !!  \frac{s-t^n}{\Delta t}K_1 +
! !!  \frac{t^{n+1}-s}{\Delta t}K_2.
! !! \f}
! subroutine rk3expo_step(ode,unew,t,unow,step_diag)
!  real(wp),     intent(in)    :: t
!  class(c_stv), intent(in)    :: unow
!  class(c_stv), intent(inout) :: unew
!  class(c_ode), intent(inout) :: ode
!  type(t_ti_step_diag), optional, intent(inout) :: step_diag
! end subroutine rk3expo_step

!-----------------------------------------------------------------------

end module mod_spexpo

